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^^ Abstract 

f) • We introduce a simple method to estimate the system parameters in continuous 

dynamical systems from the time series. In this method, we construct a modified sys- 
tem by introducing some constants (controlling constants) into the given (original) 
system. Then the system parameters and the controlling constants are determined 
by solving a set of nonlinear simultaneous algebraic equations obtained from the 
relation connecting original and modified systems. Finally, the method is extended 
to find the form of the evolution equation of the system itself. The major advantage 

^Sj , of the method is that it needs only a minimal number of time series data and is ap- 

^D I plicable to dynamical systems of any dimension. The method also works extremely 

well even in the presence of noise in the time series. This method is illustrated for 

t::;j- . the case of Lorenz system. 
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j^ ■ Time series analysis (both vector and scalar) is of considerable relevance [1,2,3] 

to physical, chemical and biological systems as they very often exhibit tem- 
poral chaotic motions. The main objective of time series analysis is to study 
the detailed structure of the evolution equation of the underlying dynamical 
system. This includes the number of independent variables, the form of the 
flow functions and parameters involved in the system. In this Letter, the study 
is focussed on the last aspect, namely, estimating the system parameters from 
the time series when partial information about the system is known (the num- 
ber of independent variables and the form of the flow functions), and then 
the study is extended to predict the form of the flow function itself when it 
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is not known. A few methods have been proposed [1,2,3] in the hterature for 
predicting the form of the flow functions. These include two point boundary 
value problem approach [4], Euler integration approach to odes [1] and modi- 
fications [5]. Also recent literature shows that the methods already proposed 
for estimating the system parameters are based on the concept of synchro- 
nization [6,7,8,9,10], Bayesian approach [11,12], least squares approach [13] 
and by successive integration method [14]. In this Letter, a very simple and 
efficient method for estimating the system parameters as well as the form of 
the flow functions of continuous dynamical systems from the vector time series 
is developed using the concept of controlling chaos [15,16,17], which can also 
in principle be extended to scalar time series. In our method we construct a 
modified system by inclusion of certain constants (controlling constants) in 
the given original system so that the evolution of the modified system is con- 
trolled to an equilibrium point. Then we find the dynamical relation between 
the original and modified systems and thereby determine the unknown system 
parameters and the controlling constants. After accomplishing this task, the 
method is extended to determine the form of the evolution equations (fiow 
functions) itself for the system from which the time series was collected. This 
method is applicable to time series obtained from a continuous system of any 
dimensions and is also well suited for discrete dynamical systems as shown in 
ref. [18] earlier. The method can also be used for the time series which contains 
considerable amount of noise. Further this method can be used in the field of 
controlling chaos to find the exact values of controlling constants to make the 
given chaotically evolving system to be controlled to a required equilibrium 
point. 

Consider an arbitrary A^-dimensional continuous chaotic dynamical system 
(the original system). 



Xi = fi{Xi,X2,...,XN;p), (1) 

where i = 1,2, 3, ...A^, and p denotes the system parameters of dimension M to 
be determined. Here we assume that the function / is sufficiently smooth. Let 
us construct a modified continuous dynamical system (the modified system) 
as 



yi = fi{yi,y2,---,yN;p) + f^i, i = i,2,...,N, (2) 

where Kj's are constants (controlling constants). The crucial idea here is that 
the Jacobian matrix which determines the stability of the equilibrium point is 
the same for both the cases, namely the original and the modified systems. In 
fact, the inclusion of Kj's in Eq. (2) makes the modified system to have an equi- 
librium point (either stable or unstable) which is effectively different from the 
equilibrium point (unstable) of the original system eventhough the Jacobian 



matrix is not affected as stated above. From an examination of many maps 
and ffows we have found tliat tliere is in general a possibility of making the 
modified system (2) to exhibit a stable equilibrium point by suitable inclusion 
of constants /tj eventhough the original system (1) evolves chaotically. 

Let Ui{t) be the deviation, that is, Ui(t) = yiit) — Xi(t), of the trajectory 
of the modified system from that of the original system due to the effect of 
controlling constants (kj's). After substituting the above relation in Eq. (2) 
and making use of Taylor expansion, we obtain 



^f. ^f^ 



1 JX, JX, 92/, 



+ oJ2J2 ^J^k 



^ j=lk=l OXjOXk 



+ ■■■, 2 = 1, 2,..., AT. (3) 



X 



Here the noteworthy point is that the above equation contains x whose evo- 
lution is characterized by Eq. (1). Thus, while obtaining the solution to Eq. 
(3) one has to solve Eqs. (1) and (3) simultaneously. 

Let us now consider the time evolution of the variables Xi{t), yi{t) and 'Uj(t), 
statisfying Eqs. (l)-(3), respectively, at discrete time intervals. For this pur- 
pose, we introduce the notation x^ = Xi{jAt), yl = yi{jAt), and u} = 
Ui{jAt), where At is a small time interval and j = 0,1,2,3,.... With this 
notation, the relation between the original and modified systems at definite 
intervals can be written as 



yO)=^<J)^^)^ J =0,1,2,... (4) 

After the transient time kAt, let the modified system approaches an equilib- 
rium point y* and hence the above equation becomes. 



xp) + Mp) =y*, j = k + l,k + 2,k + 3, ... (5) 



Let 4 ,4 , ..., 4™ be the m set of data points of the given chaotic time se- 
ries sampled at the time interval At from the original system (1). Substituting 
this in the above equation, we get 



z^ + m1") = y*, n = 0, 1, 2, ..., (m - 1). (6) 



In the above relation, it is instructive to note that we need not bother about 
the transients because the time series data is collected only after sufficient 
transient time. So, the discrete time index (n) now can start from 0, that is 



from the first data of the time series. It may be noted that here, -Uj's are func- 
tions of the parameters p and controlhng constants /tj, since the derivatives 
of -Uj's are having specific functional relations with the same parameters and 
controlling constants through Eq. (3). Also, Mj s are obtained by solving Eq. 
(3) numerically with time interval At (or submultiples of At). In general, for a 
given set of time series data collected from the system (1) for a particlular set 
of parameters, the right hand side of Eq. (6), that is, the equilibrium point, 
depends on the value of controlling constants (kj's) which can be varied by 
means of m^^ . This freedom allows us to make any desired point in the region of 
the attractor of the system as equilibrium point (y*) by choosing the values for 
u\ accordingly. For example, one can easily have z^ as an equilibrium point 
by setting u\ = in Eq. (6). Similarly, an arbitrary point qi in the attractor 
can be chosen as an equilibrium point if we start with u\ = qi — zl . 

Now it is possible to obtain the required number (M + N, N for /tj's and 
M for p's) of functional relations (implicit) between the unknown parameters 
and controlling constants through Eq. (6). Once we have the required number 
of functional relations, the next task is to solve them for the unknowns, that 
is for Kj's and p's, using a suitable numerical technique, for example by the 
globally convergent Newton's method [19]. 

After estimating the values of the parameters, its accuracy can be checked as 
follows. Compare the equilibrium point (for example y* = zi if we assume 
M^ = in Eq. (6)) assumed to be exhibited by the modified system in the 
above procedure with the equilibrium point calculated from Eq. (2) with the 
values of parameters determined by the above method. The degree of close- 
ness of these equilibrium points gives a measure of accuracy of the estimated 
parameters. 

As an illustration to our method, let us consider the Lorenz system 



Xi = (T(a;2 — xi), (7a) 

±2 = pXi — X2 — X1X3, (7b) 

X3 = xiX2-bxs, (7c) 

where a, p and b are the unknown system (control) parameters. Then the 
modified Lorenz system can be constructed as 

yi = (^iy2-yi) + K,i, (8a) 

m = pyi - ^2 - 1/11/3 + «2, (8b) 

y3 = yiy2-bys + Ks, (8c) 

where Ki, K2 and /ta are constants to be determined which make (8) to exhibit 
equilibrium point solution while the original system (with ki = K2 = n^ = 0) 





Table 1. Convergence 


of cr, p, 6, Ki, 


K2and K3 in the Lorenz system 
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f^2 


/«3 





1.00000000 


1.00000000 


1.00000000 


1.00000000 


1.00000000 


1.00000000 


1 


25.06611000 


-53.43927900 


5.90702312 


-1.09685478 


-1.64477291 


5.66012411 
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33.88597350 


5.87195581 


4.74703140 


-1.70082489 


-2.33269811 


7.14465601 
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-24.31336050 


20.11146620 


3.81209093 


-3.00648132 


-3.72360435 


9.50270473 
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31.73494710 


37.21100300 


1.76763806 


-7.61414582 


-9.53379493 


19.84796380 
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-11.56754590 


20.78649710 


3.33689508 


-15.06851440 


-19.96489090 
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Fig. 1. Convergence of the parameters o", p and h towards their exact values in the 
Lorenz system. 

exhibits chaotic solution. Then the equation (3) for the deviations, Ui = Hi—Xi, 
i = 1,2, 3, becomes 



ill = cr{u2 — ui) + Ki, 

U2 = {p- X3)Ui - U2- {Xi + ■Ui)-U3 + K2, 
U3 = {X2 + U2)Ui + X1U2 - bus + /«3- 



(9a) 
(9b) 
(9c) 



It may be noted that the presence of Xi, X2 and X3 in (9) indicates that one 
has to solve Eqs. (7) & (9) simultaneously. 
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Let (z'^, zli, z'^), {zl, Z2, z^), ..., (zi" ^' , Z2"' ''' , z^" ''') be the time series data 



(m-l)^ 
3 y 

obtained from the Lorenz system at some arbitrary time interval with the 
sampling rate At for an unknown set of system parameters. After assuming zl 
as the equilibrium point for the modified Lorenz system (that is by assigning 
u\ = in Eq. (6), so that y* = zl ) and substituting three data points 4 , 
zl and zl (one can in fact take any three successive data points but the 
first data point has to be used as initial condition for Xi) in Eq. (6), we get 



z« + 4'^ 



.(2) 



uf^ 
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(10) 



Note that the time derivatives of UiS have specific functional relations with 
the unknown parameters {a, p and 6) and controlling constants {ki, K2 and 
^3) through Eq. (9), and hence Mj's are also functions of the same parameters 
and controlling constants. Now we have six (implicit) functional relations for 
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Fig. 2. The distribution of the values of the system parameters estimated from the 
time series which contains random dynamical noise in the Lorenz system. 



the six unknowns namely a, p, b, Hi, K2 and K3 through Eq. (10), in terms 



.(0) 



(1) 



of the equilibrium point (here z\ ) and the two sets of data points z\ and 
1, 2, 3. To solve for these unknowns, we have to obtain values of 



u 



.(2) 
i ) 
(1) 



U 



(1) ^,(1) ^,(2) ^,(2) 



.(2) 



1x3 , u\ \ U2 and M-j by solving Eqs. (7)&(9) simultaneously so 
that Eq. (10) is satisfied. This can be done, for example, by the globally 
convergent Newton's method [19] provided an initial guess is given to all the 
unknown parameters. In this example, while obtaining the value of u\ initial 
conditions should be taken as x^ = zl and u\ =0, respectively. Further, 
the evaluation of u\ needs to reset xi = z] in the numerical algorithm used 
above. Similar procedure has to be followed if any other set of successive data 



is used in place of z\ \ z\ and z\ in Eq. (10). For illustration purpose, we 
have used the numerically generated time series of the Lorenz system for the 
system parameters a = 10.0, p = 28.0 and 6 = 8/3 and solved the Eq. (10) 
by globally convergent Newton's method with an initial guess 1.0 to all the 
unknowns Ki, k,2, K3, (J, p and b. The convergence of the system parameters 
and controlling constants (k's) is shown in the Table I, which shows that 
the estimated values are indeed the exact values of the parameters at which 
the time series data of the Lorenz system is generated. We also note that the 
convergence is very rapid, which is further confirmed in Fig. (1). Note that if we 
take a different set of data points from the time series as the equilibirium point, 
the corresponding controlling contstants (kj's) will be different eventhough the 
system parameters are unaltered. 



In order to show that our method is robust, a time series which contains 
random dynamical noise of strength 10~^ was additionaly introduced in the 
above Lorenz model. In this case, the distribution of the values of the system 
parameters estimated using 1000 data points is shown in fig. (2). Also in the 
case of random observational noise of same strength, the distribution of the 
values of the parameters is shown in fig. (3). The peaks about the true values 
of the parameters in the above two figures indicate that our method works 
extremely well even in the presence of dynamical noise or in the presence 
observational noise in the time series. 
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Fig. 3. The distribution of the values of the system parameters estimated from the 
time series which contains random observational noise in the Lorenz system. 

We have also carried out a similar analysis for the autonomous Chua's Circuit 
[16] in its dimensionless form, 



Xi = a{x2 — Xi — h{xi)), 
a^2 = a;i - X2 + X3, 

X3 = -(3X2, 

where 

h{x) =bx + 0.5(a - b){\x + 1| - |a; - 1|), 

a and f3 are parameters to be determined, and for the Rossler system. 



(11a) 
lib) 
die) 



X^ = -{X2 +X3), 

a;2 = xi + ax2, 

X3 = b+ {xi -c)x3. 



(12a) 
(12b) 
(12c) 



where a, b and c are the parameters to be determined. In both the cases, we 
are able to obtain the correct values of the parameters as in the case of the 
Lorenz system. 

Next, we wish to point out that it is also possible to extend the analysis to 
predict the flow function of the system itself in principle, by assuming a poly- 
nomial form (here as an illustration) for the functions fi{xi,X2, ...Xn', p) in the 
right hand side of Eq. (1), with unknown coefficients, and solving sufficient 
number of Eqs. (6) to determine them. One can as well choose other types 
basis functions and our procedure is equally applicable here also. To illustrate 
this idea, let us consider the same time series data used in the above exam- 
ple of Lorenz system and assume that the form of the underlying dynamical 
equations is unknown. Let us assume a general quadratic form for the flow 
function (dynamical equations) as 



xi = ciXi + C2X2 + 03X3 + C4X1X2 + C5X2X3 + ceXiXs, (13a) 

2^2 = c^Xl + C8a;2 + C9X3 + Cioa;ia;2 + Ciia;2a:3 + £12^1X3, (13b) 

xz = C13X1 + C14X2 + C15X3 + 016X1X2 + C17X2X3 + C18X1X3, (13c) 

where q's, i = 1,2, ...18, are parameters to be determined from the time series 
data. Again following the method suggested above, we write the form of the 
modified system as 

yi = civi + C2I/2 + C3I/3 + C4I/1I/2 + C5I/2I/3 + C6I/1I/3 + Ki, (14a) 

?/2 = CjVi + C8t/2 + C9I/3 + CiQyiy2 + Ciit/2l/3 + C12I/1I/3 + /«2, (14b) 

?/3 = Cl3|/1 + Ci4|/2 + Ci5|/3 + C16I/1I/2 + C17I/2I/3 + C18I/1I/3 + /«3, (14c) 

where ki, k,2 and ^ are again the controlling constants to be determined so 
that the above modified dynamical system exhibits a stable equilibrium point. 
For this system, Eq. (3) becomes 



Ul = Ki+ CiUi + 02^2 + C3U3 + (C4X2 + C6X3 + C4U2)Ui 
+ (C4X1 + C5X3 + C5M3)m2 + (C5X2 + CeXi + CqUi)u3, 

U2 = K2 + C-jUi + C8-U2 + £9^3 + (CioX2 + C12X3 + CiqU2)ui 

-H'CioXi + C11X3 + Cii-U3)-U2 + (C11X2 + C12X1 + Ci2-Ui)m 

+ C13M1 + C14M2 + C15M3 + (ci6a:2 + Ci8a:3 + 016^2)^1 



+ 

«3 = /«3 



+ (Ci6Xi + C17X3 + CnUs)u2 + (C17X2 + C18X1 + Ci8Mi)m3. (15c) 

In order to determine the values of the parameters and controlling constants 
(cj, i = 1,2, ...18, Ki, K2 and k^), we have to solve 21 algebraic equations which 
are actually constructed by making use of eight successive time series data in 
Eq. (6) with assumption that the first one is the equilibrium point. Then the 
required equations will be 

4^')+Mp)=zf, ^ = 1,2,3 and J = 1,2,.. .,7. (16) 



Again we follow exactly the steps mentioned earlier for the Lorenz system 
and obtain the values of the unknown parameters. The distribution of the 
values of the parameters (q's, i = 1,2, ...,18) obtained by solving Eq. (16) 
using 3000 data of time series is shown is Fig. (4). And the values are found 
to be distributed around {ci}\^ = {-10, 10, 0, 0, 0, 0, 28, -1, 0, 0, 0, -1, 
0, 0, —2.67, 1, 0, 0}. Substituting these values in Eq. (13) we obtain the flow 
function (evolution equation) of the form 

xi = CiXi + C2X2, (17a) 

X2 = C7X1 + C8X2 + C12X1X3, (17b) 

a:3 = ci5a;3 + ci6a;iX2, (17c) 
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Fig. 4. The distribution of the values of the system parameters of the flow functions 
given by Eq. (13), estimated from the time series of Lorenz system 
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with a = — ci = C2 = 10.0, p = cj = 28.0, h = — C15 = 2.67 and the remaining 
constants Cig = 1.0 and Cg = C12 = —1, which is nothing but the original 
Lorenz system. 

Finally, we would like to point out that the procedure outlined above also gives 
a method to obtain the values of the controlling constants (kj) for a chaotic 
system to be controlled to a desired equilibrium point if the form of evolution 
equation is known. 

To conclude, we have developed a very simple as well as useful method for 
estimating the unknown system parameters of the continuous dynamical sys- 
tems of any dimensions from the vector time series, if partial information is 
known, namely the form of the dynamical equations . It has also been ex- 
tended to obtain the form of the system equation itself at least in the case of 
polynomial forms. Both the methods have been illustrated by means of vector 
time series collected from the Lorenz system, while further confirmation has 
been made with other systems including the Chua's Circuit and Rossler sys- 
tems. We have also checked that the method is robust againt dynamical and 
observational noise and that the procedure exhibits a rapid convergence. 
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